import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.stats import norm

# 解决图形显示中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 生成正态分布样本
a = np.random.normal(loc=0, scale=1.0, size=1000000)
print('正态分布样本:', a)

# 样本均值的无偏估计
b = np.mean(a)
print('样本均值的无偏估计：', b)

# 样本方差的无偏估计
c = np.std(a, ddof=1)
print('样本标准差的无偏估计：', c)

# 绘制正态分布图
plt.hist(a, density=True, bins=30)  # 归一化直方图（用出现频率代替次数），将划分区间变为 30（默认 10）
x = np.linspace(-3, 3, 50)  # 在(-3,3)之间返回均匀间隔的50个数字。
plt.plot(x, norm.pdf(x), 'r-')
plt.title('正态分布图')
plt.show()


# 正态分布均值区间估计
def confidence_interval_u(data, sigma=-1, alpha=0.05, side_both=True):
    xb = np.mean(data)  # 均值
    # s = np.std(data, ddof=1)  # 标准差
    # sigma已知，服从标准正态分布
    z = stats.norm(loc=0, scale=1.0)
    if side_both:  # 求双侧置信区间
        tmp = sigma / np.sqrt(len(data)) * z.ppf(1 - alpha / 2)
        m = xb - tmp
        n = xb + tmp
        return m, n
    else:  # 单侧置信下限或单侧置信上限
        tmp = sigma / np.sqrt(len(data)) * z.ppf(1 - alpha)
        m = xb - tmp
        n = xb + tmp
        return {'bottom_limit': m, 'top_limit': n}


def confidence_interval_sigma(data, mu=-1, alpha=0.05, side_both=True):
    # xb = np.mean(data)
    # s_square = np.var(data, ddof=1)
    sum_tmp = 0.0
    for i in data:
        sum_tmp = sum_tmp + (i - mu) ** 2
    if side_both:  # 双侧置信区间
        m = sum_tmp / stats.chi2.ppf(1 - alpha / 2, df=len(data))
        n = sum_tmp / stats.chi2.ppf(alpha / 2, df=len(data))
        return m ** 0.5, n ** 0.5
    else:  # 单侧置信下限或单侧置信上限
        return {'bottom_limit': sum_tmp / stats.chi2.ppf(1 - alpha, df=len(data)),
                'top_limit': sum_tmp / stats.chi2.ppf(alpha, df=len(data))}


h = confidence_interval_u(a, 1)
k = confidence_interval_sigma(a, 0)
print('样本均值的区间估计:', h)
print('样本标准差的区间估计:', k)